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We present a comparative study of the rotational characteristics of various 
molecule-doped 4 He clusters using quantum Monte Carlo techniques. The the- 
oretical conclusions obtained from both zero and finite temperature Monte Carlo 
studies confirm the presence of two different dynamical regimes that correlate with 
the magnitude of the rotational constant of the molecule, i.e., fast or slow rotors. 
For a slow rotor, the effective rotational constant for the molecule inside the he- 
lium droplet can be determined by a microscopic two-fluid model in which helium 
densities computed by path integral Monte Carlo are used as input, as well as by 
direct computation of excited energy levels. For a faster rotor, the conditions for 
application of the two-fluid model for dynamical analysis are usually not fulfilled 
and the direct determination of excitation energies is then mandatory. Quanti- 
tative studies for three molecules are summarized, showing in each case excellent 
agreement with experimental results. 

1 Introduction 



Droplets of 4 He provide a unique quantum environment that constitutes an ultra- 
cold and gentle matrix for high resolution spectroscopy 1 and for investigating the 
nature and dynamic consequences of quantum solvation for a wide range of im- 
purities doped into superffuids. 2 One of the most unusual features deriving from 
those experiments is the apparent free rotation of small molecules embedded inside 
these bosonic clusters at temperatures T ~ 0.4 K. 3 Rotational spectra obtained 
for a series of molecules possessing gas phase rotational constant B in the range 
0.01 — 50 cm , ' appear to fall into two dynamical regimes. The heavier molecules 
(slow rotors) show a reduction in rotation constant of ~ 60 — 80%, while lighter 
molecules (fast rotors) show a much smaller reduction of — 30%. In this article, 
we focus on the behavior of the three molecules SFg, OCS, and HCN. The first two 
belong to the first dynamical regime of slow rotors, showing a similar reduction of 
B/Bq ~ 37% in both cases. 1 The third molecule, HCN, is a much faster rotor and 
the experimentally observed reduction of rotational constant B/Bq is only 19%. 4 
For each of these molecules, accurate pair interaction potentials with helium ex- 
ist, rendering quantitative calculations meaningful and reliable tests of quantum 
theories for molecular rotation in helium droplets. 
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2 Methodologies 

Our theoretical studies of molecule-doped helium clusters are based both onT = 
methods diffusion Monte Carlo (DMC) and the projector Monte Carlo extension of 
this for excited states (POITSE), and on the finite temperature path integral Monte 
Carlo (PIMC) approach with incorporation of full quantum exchange symmetry. 3 
With DMC based techniques we evaluate directly the ground and excited energy 
levels of the doped cluster, focusing in this study on rotational excited states. 
With the PIMC calculations we are able to compute finite temperature helium 
densities. These can then used as input for a two-fluid dynamical model that 
combines notions of adiabatic following with hydrodynamic estimates of superfluid 
response to molecular rotation. 3 We provide in this section a short summary of how 
the quantum Monte Carlo methods arc implemented for study of systems composed 
of a single molecule treated as a rigid body, together with N helium atoms. 



2.1 Diffusion Monte Carlo 

The diffusion Monte Carlo technique for solving the Schrodinger equation for many- 
body systems is based on the imaginary time (r = it/h) Schrodinger equation 

^ffi = £ DjV^m [V(H) E ref ] (1) 

3 

where 1Z is a vector in the A/"-dimensional space, Dj — h 2 /2mj if the j th degree of 
freedom corresponds to a translation, and Dj = Bj = h 2 /2Ij if it corresponds to a 
rotation. In the above equation, E re f is a constant defining the zero of the absolute 
energy scale, and V is the potential. This formulation implies the use of Cartesian 
coordinates for the atom-like particles and for the center of mass of the rigid body, 
and the use of rotational angles around the principal axes of the rigid body. This 
diffusion-like equation is solved using random walks and a short time approximation 
for the Green's function. 5 The effects of statistical noise on the results can be 
drastically reduced by introducing a guiding (or "importance sampling") function 
that approximates the true solution of Eq.(l), and by then rewriting the 
diffusion- like equation for the product function f(TZ) — ^ (IZ)^ t {TZ) : 

^ffi = E i D ^m) D j V j [fWFjill)}} mil) E ref ] f(R). (2) 

3 

Here Ei(R) = *t(^) _1 ^*t(7?.) is the local energy, and F (TZ) = Vj In |*t(^)| 2 
is the quantum force that controls the drift terms (second set of terms on right- 
hand side of Eq. (2)). We have recently shown 6 that one can use this importance 
sampling scheme for all degrees of freedom, including the rotation of the molecule 
when treated as a rigid body. The average of the potential (in unbiased walks) or of 
the local energy (in biased, or "importance sampled" walks) over the random walk 
yields the ground state energy of the system. Whenever a short-time expansion 
for the Green's function is used, it is necessary to check that any time step bias is 
removed, or is at least smaller than the statistical noise. 



Copyright © 2002 by World Scientific, reprinted with permission 2 



With the incorporation of importance sampling into DMC it becomes possible to 
compute the helium density, or at least a mixed density (^t\p\^), by simply binning 
the position of the helium atoms in a frame attached to the molecule during the 
random walk. It is also possible to analyze the extent of adiabatic following of the 
molecular motion by the helium density, by artificially suppressing the rotational 
kinetic term of the dopant molecule in the Hamiltonian operator. 3,7 ' 8 Complete 
adiabatic following implies that the helium density in a frame attached to the 
molecule must be identical in both cases, i.e., with and without the incorporation of 
molecular rotation. We have also proposed 3 ' 8 a straightforward method to quantify 
the extent of the adiabatic following. This can be done by evaluating the ratio 
of the densities along privileged axes in the molecule-helium system, with and 
without rotation: Q(r) = ([p sa ddie/pmin] no rot ) / ([psaddie//Omin] rot )- When adiabatic 
following is significant, Q(r) is close to one. In contrast, a value of Q(r) close to 
zero indicates weak or minimal adiabatic following. 

Just as in analysis of experimental spectra, 1 the rotational constants can be 
extracted by fitting the energy differences between different rotational states. For 
spherical and linear rotors, the first energy difference between ground and first ex- 
cited states is equal to IB. This approach requires the explicit computation of 
excited state energies. Such computations can be done in an approximate way by 
making use of the widely employed fixed node approximation. In this approxima- 
tion, one imposes a predefined nodal surface on the excited state wave function by 
use of a trial wave function bias, Eq. (2). However, unless the trial function node 
is determined by symmetry considerations, great care must be taken when using 
such a method. As we have seen for the dimer HCN-He, 6 ' 9 a non-physical node can 
easily lead to spurious energy levels. Exact calculations for N = 1 "dimers" have 
been useful for calibrating fixed node calculations. A more reliable approach to cal- 
culation of excited states is provided by the intrinsically more powerful and exact 
method of Projection Operator Imaginary Time Spectral Evolution 10 (POITSE). 
This method allows the computation of excitation energies without imposing any 
fixed node constraints. The POITSE approach makes use of the inverse Laplace 
transform of a "projected" correlation function that leads to the spectral function 

k(E) = J2 \(HM^f)\ 2 S(E - E f + E). (3) 
/ 

Different excited energies Ef are accessed by choosing the projector A such that 
the overlap {tpo\A\ip f) is non- vanishing for one or a few excited states ipf. The 
evaluation of the correlation decay that constitutes the transform of Eq. (3) is done 
by making DMC side walks on the Metropolis sampling of a ground state density. 10 
The numerical implementation used in the studies we summarize here has been 
described elsewhere. 9 ' 11 It is based on performing the DMC side walks with pure 
branching, rather than the pure weights version that was used originally. 10 This 
allows for the computation of excited states of much larger systems than studied 
previously. 

For both fixed node and POITSE calculations, one needs to have some infor- 
mation about the nodal structure of the excited state of interest. In the case of 
fixed node, this is obviously extremely important since as mentioned above, the 
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results are completely dependent on the choice of this nodal surface. 9 In contrast, 
with POITSE, one can reach the exact energy level even with a "wrong" projector. 
This is nicely illustrated by the different behavior in 4 He at of the three molecules 
studied here. Examination of the symmetry of rotational spectra obtained in he- 
lium droplets (spherical top for SF 6 and linear molecule for OCS), as well as the 
theoretical study of rotational excitations, made for the reduced dimensionality 
"diatomic" SF 6 model 12 indicates that the relevant excited energy levels are those 
in which the rotational energy is located primarily on the molecule, i.e., J ~ j. 
One corresponding wave function for this situation in the case of a free molecule is 
the Wigner function \j = 1, k = 0, m = 0). One can use this Wigner function as a 
trial nodal structure or, alternatively, as a projector for POITSE. In the fixed node 
approximation, the results are valid only insofar as this particular nodal constraint 
is valid, i.e., physically correct. In contrast, application of the POITSE method 
with this Wigner projector leads to the correct energy levels no matter how weakly 
the true nodal structure in the cluster resembles that of the free molecule. 

For small N, these rotational constants derived from direct computations of the 
rotational energy levels of the doped cluster can be usefully compared with values 
obtained from rigid coupling approximations for the total helium density. This 
extreme dynamical approximation can also be described as a form of adiabatic 
following, although it is important to realize that the consequences of following 
with rigid coupling of the total helium density arc different from both the adiabatic 
following of a hydrodynamically responding fluid, 3 ' 13 and the adiabatic following 
analysis for a two-fluid decomposition of the density. 3 ' 8 

2.2 PIMC, two- fluid, and hydrodynamic models 

The path integral Monte Carlo approach allows direct calculation of thermal aver- 
ages of observables, with incorporation of the boson permutation symmetry for He 
and without introduction of any trial function bias. It is currently also the only 
numerical method capable of directly addressing superfluidity at non-zero temper- 
atures. 

In PIMC one computes the average of a quantum operator over the thermal 
density matrix: 



Here p(TZ, TZ'; (3) = (TZ'\e^ l3H \1Z) is the thermal density matrix in the position rep- 
resentation, Z = J dlZp(JZ, 1Z; (3) is the partition function, and (3 — l/k^T. The 
path integral Monte Carlo approach to evaluating Eq. (4) starts from the discrete 
Feynman path-integral expansion of p(7l, TV \ 0) in terms of high temperature com- 
ponents p(7Zi,7Z2;t), where r = /3/M constitutes an imaginary time step. Use 
of the discrete path integral expansion relics on the ability to find accurate high- 
temperature density matrices p{TZ\, TZ2] t) that render evaluation of the consequent 
highly multi-dimensional integrals worthwhile, given the intractatibility of finding 
representations for the low temperature density matrix of an interacting quantum 
system, p(R, TV;/3). For 4 He systems, the density matrix must also be symmetrized 
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with respect to particle permutation exchanges, increasing the computational com- 
plexity. This can be done with multi-level Metropolis sampling schemes that pro- 
vides simultaneous sampling of the permutation and configuration space, allowing 
direct numerical study of superfluid properties at finite temperature, in addition to 
structural and energetic properties. 14 

Path integral calculations provide estimators for thermal expectation values 
(Eq. (4)). Linear response estimators for global superfluid fractions can be defined 
for both homogeneous and inhomogeneous systems, bulk or finite in extent. 14 Thus 
quantification of a superfluid on a nanoscale is accessible from the global superfluid 
fraction for a pure helium droplet. The study of molecules in helium clusters has 
raised the question of what the local superfluid density around a molecule is, and 
how this is distributed over the quantum solvation structure induced by the molecu- 
lar interaction. 3 Recent work in our group has proposed a local superfluid estimator 
in terms of the length of the Feynman paths. 15 This has been used to establish 
a local two-fluid model of the quantum solvation structure, based on the finding 
from PIMC calculations that within the first solvation shell, a local nonsupcrfluid 
density is induced by the molecular interaction with helium. 3 ' 15 The extent of this 
molecular-interaction induced nonsupcrfluid density depends on the strength of the 
interaction potential. 

This local two-fluid model has provided an atomic scale analysis of superfluid 
solvation, 15 and has also been used to analyze the rotational dynamics of the 
molecule. 3 The key additional concept for analysis of molecular rotation is the 
notion of adiabatic following of some or all of the local helium density with the 
molecular rotation. This adiabatic following, already encountered and quantified 
in the T = diffusion Monte Carlo calculations discussed above, leads to very dif- 
ferent consequences for the local nonsuperfluid density and for the local superfluid 
density. For the local nonsupcrfluid density, adiabatic following implies an effective 
rigid coupling to the molecular rotation, while adiabatic following of the local su- 
perfluid density can give rise to hydrodynamic added inertia, 3,5 analogous to the 
well-known added mass of macroscopic hydrodynamics. As discussed in detail in 
Ref.3, this results in a microscopic analog of the Andronikashvili experiment, with 
the additional distinction that the molecular rotation is quantized and does not 
provide a continuously varying driving force, unlike a macroscopic rotating probe. 
This results in additional angular momentum constraints when the two-fluid analy- 
sis of molecular rotation is applied to interpretation of spectroscopic experiments. 3 
An alternative pure hydrodynamic model is obtained as the limiting case of the 
two-fluid model when the entire density is assumed to be superfluid and undergoes 
adiabatic following. 3 ' 13 

3 Applications 

In this section we present a summary of the various quantum Monte Carlo appli- 
cations made to date for molecules embedded in 4 Hc clusters. We focus on three 
different molecules: SF 6 , a slow spherical rotor, OCS, a slow linear rotor, and HCN, 
a fast linear rotor. Table I summarizes the corresponding rotational constants in 
the gas phase (Bo) and in helium clusters (B), and also gives the molecular masses 
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used in the calculations. 



Table 1. Experimentally rotational constants in gas phase Bo and inside helium droplets B, in 
cm -1 . Also reported are the corresponding masses, in amu. 



Molecule 


So 


B 


Rcf. 


mass 


SF 6 


0.091 


0.034 


16,17 


146.0288 


ocs 


0.20 


0.073 


18 


60.07455 


HCN 


1.47 


1.20 


4 


27.01104 



3.1 SF 6 

The first experimental rotational spectrum was obtained for the spherical top 
SF6. 16 ' 17 For this slow rotor, the reduction of the B value inside a helium cluster is 
relatively large, with B/Bo ~ 36%. This molecule was also the first for which the- 
oretical studies of the rotational dynamics were made. 3 ' 7,12 Using the anisotropic 
potential energy surface of Pack and co-workers, 19 SF 6 He at has been studied both 
with PIMC, 3 and with DMC in the fixed node approximation. 7 A POITSE study 
of a diatomic model of SF 6 has also been made. 12 

The nonsupcrfluid density of helium computed via PIMC is located in the 
first solvation shell and shows angular modulation structure resulting from the 
anisotropic interaction of helium with the octahedral SFg molecule. 2 ' 3,15 The ener- 
getic criterion for adiabatic following is fulfilled in this case, and application of the 
two- fluid dynamical model leads to a value of B = 0.033 cm^ 1 in very good agree- 
ment with the experimental value. Hydrodynamic calculations employing PIMC 
densities calculated at the experimental temperature of T = 0.3 K show negligible 
added mass contribution from either the superfluid or total density. 3 The average 
of the Q-factor for the dimer SF 6 -He is found to be around 0.7. 8 This value is large 
but not equal to unity, indicating that the adiabatic following of the helium density 
with the molecular rotation is substantial but not complete in this system. The 
Q-factor is similarly less than unity in larger clusters, 3 implying that there is only 
partial adiabatic following by the density within the first solvation shell. This is 
consistent with the detailed dynamical analysis of the two-fluid model. 3 

The extraction of the rotational constant B has also been performed by direct 
computation of energy levels. J = 1,2,3 levels for clusters with N = 1 to 20 
helium atoms were calculated in fixed node, employing the free molecule nodal 
surfaces. For all sizes, the levels can be beautifully fit to the spherical top formula 
BJ(J + 1). The resulting B value decreases monotonically from the gas phase 
value as N increases from unity until it reaches a plateau at N — 8 (Fig. 1). An 
asymptotic value of B = 0.035 cm -1 is obtained, in very good agreement with the 
experimental value extracted from spectra in large 4 He droplets (Table 1). 

In the POITSE study for SF@, the octahedral SFq molecule was represented 
by an effective diatomic. 12 The use of the original numerical implementation of 
POITSE employing continuous weights in the DMC sidewalks limited this study to 
N = 1,2 and 3 helium atoms. Nevertheless, the extracted B values are similar to 
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the values obtained from the full dimensionality DMC calculations within the fixed 
node approximation (see discussion in Ref. 3), indicating that the free molecule 
nodal structure used is indeed physical in this particular case. 

1.2 
1.0 
0.8 
« 0.6 
0.4 
0.2 
0.0 

1 2 3 4 6 8 10 12 15 20 25 

N 

Figure 1. Evolution of B/Bq with the number of helium atoms N for SFg (circles) and HCN 
(triangles). The experimental values in large droplets, 36% (SFg) and 81% (HCN), are indicated 
by the arrows on the right hand axis. For SFg the theoretical results were obtained within the 
fixed node approximation, 7 and for HCN from POITSE calculations. 9 



3.2 OCS 

The second dopant molecule studied by quantum Monte Carlo methods is OCS, 
which is a linear rotor possessing a gas phase rotational constant Bq twice as big as 
that of SF 6 . The increase of the moment of inertia measured in helium droplets 18 
leads to a ratio B/B ~ 37 % that is nevertheless similar to the ratio for SF 6 . 
Because the anisotropy of the OCS-He potential 20 is very pronounced, we expect 
that adiabatic following will also be strong for this system, even though the rotor 
is somewhat faster than SFq. Path integral computation of the local supcrfluid 
and nonsuperfluid solvation densities show that the energetic criterion for adiabatic 
following is indeed fulfilled. 3 Computation of the rotational constant using the two- 
fluid dynamical model leads to a value B = 0.067cm _1 when angular momentum 
constraints are applied, 3 in good agreement with the experimental value reported 
in Table 1. For this system, the Q-factor computed by DMC is also around 0.7, 
from which we can conclude that for OCS, adiabatic following of the helium density 
is significant but not complete. 

When a quasi-adiabatic rigid coupling analysis 3, 21,22 is made for a small num- 
ber of helium atoms (N < 7) , this shows that as for SFg , the rotational B constant 
decreases monotonically with N, and that a saturation will probably occur for 5 or 
6 helium atoms. However, as pointed out in Ref. 3, and confirmed in Ref. 23, this 
saturation cannot be seen within the rigid coupling approximation, which continues 
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to show a monotonic decrease of B as the cluster size increases. In contrast, the 
POITSE approach allows calculation of the rotational excitation without any dy- 
namical approximation. POITSE calculations have recently allowed the saturation 
at small N values to be demonstrated explicitly. 22 A projector based on the free 
molecule Wigner function, \j = 1, k = 0, m = 0) produces a POITSE decay which 
contains only a single exponential term. This is easily Laplace inverted using the 
maximum entropy approach. The position of the single peak obtained, which we 
equate with 2B, moves towards a smaller value as N increases from 1 to 6. The 
peak position then remains unchanged as N is further increased to 10 and to 20. 
The extracted saturation B value of 0.07 cm -1 is in good agreement with the ex- 
perimental value measured in much larger clusters (see Table 1). This saturation 
at a small number of helium atoms corresponding to only a fraction of the first 
solvation shell is similar to the situation described above for the octahedral SFg 
molecule. 7 It similarly implies that for OCS a fraction of the first solvation shell 
density undergoes adiabatic following characterized by an effective rigid coupling 
of only that density fraction to the molecular rotation. 

3.3 HCN 

The last dopant molecule we discuss here, HCN, is also a linear rotor, but a much 
faster one than OCS. The value of Bo for HCN is more than one order of magni- 
tude larger than the values for SFg and OCS. The reduction of rotational constant 
observed experimentally for HCN is also much smaller than the corresponding re- 
duction for the two previous molecules. Thus, for HCN the ratio B/B$ is ~ 81%. 4 
Path integral analysis concluded that both the two-fluid dynamical model and hy- 
drodynamic models are not applicable here because of lack of efficient adiabatic fol- 
lowing of the HCN rotation. 3 This is not surprising, given the very weak anisotropy 
of the interaction potential for HCN with a single helium atom 24 and the relatively 
high value of the gas phase rotational constant B . 

Confronted with the impossibility of applying two-fluid or hydrodynamic models 
in this case of a fast rotor, energy level calculations have been performed to directly 
access the rotational excitations in helium clusters. 9 These have been made with 
both fixed node and with the exact POITSE method, and the effective rotational 
constant B extracted in each case by fitting energy level differences. By comparison 
with levels obtained from basis set calculations, 9 ' 24 it was found that for HCN at 
the very smallest cluster size (N — 1), the nodal structure of the free molecule is 
not physical and leads to an erroneous energy level that does not appear in the 
basis set calculations. We have proposed that this feature derives from an unusual 
coincidence between the ground state energy and the potential barrier in the HCN- 
He system. 9 

Since the nodal structure of the excited state is not easily approximated for this 
light molecule, the POITSE approach is certainly the method of choice here. We 
have shown that the spectral inversion of exponential projector correlation decays 
obtained from the same projector as that used for OCS (see above) leads now to 
spectra with two peaks. 9 Excellent agreement is found while comparing the posi- 
tions of these two peaks for HCN- He with results of basis set calculations, implying 
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that this projector now overlaps with two excitations neither of which possesses 
the free molecule nodal structure. This difference in the number of excitations 
seen in the POITSE spectra obtained from the same projector for the OCS and 
HCN molecules derives from the fact that these two systems lie in very different 
dynamical regimes, for which quite different nodal structures appear to apply. 

The size dependence of the calculated B value for HCN is compared with that 
for the heavier SF 6 molecule in Fig. 1. It is evident that molecules in the two 
different dynamical regimes show different saturation behavior. While the heavy 
molecule shows a monotonic decrease in B / Bq to its saturation value, the light 
molecule first undershoots and then gradually increases to its saturation value at 
a significantly larger value of N. This slower saturation behavior is believed also 
to be related to the complex interplay between rotation and potential hindrance in 
the HCN-He system. 9 In both cases it is nevertheless evident that saturation to 
the experimental value of B/B occurs at a cluster size N that is far smaller than 
the experimental cluster sizes (N ~ 4000 1 ), confirming that the bosonic quantum 
solvation effect on dopant rotations in 4 He is primarily a local effect. 

4 Conclusion 

We have presented a comparative study of several doped helium clusters using 
quantum Monte Carlo techniques. For the heavier two dopants presented here, 
SF 6 and OCS, the path integral calculations show that the quantum solvation 
in superfluid He exhibits a molecule-induced nonsuperfiuid density in the first 
solvation shell. The extent of this nonsuperfiuid density is directly related to the 
strength of the impurity-helium interaction potential. For HCN, the interaction 
potential is sufficiently weak that a nonsuperfiuid density cannot be consistently 
defined. With both finite and zero temperature techniques, we are able to calculate 
the extent of adiabatic following of molecular rotation by a fraction of the helium 
density. The adiabatic following calculated for the two heavier molecules, SF 6 and 
OCS, is quite similar and strong, whereas for HCN the extent of adiabatic following 
is very small. This strong difference between the heavy and light molecules has the 
consequence that neither two-fluid, pure hydrodynamic, nor quasi-adiabatic rigid 
coupling models apply for the computation of the rotational constant of HCN in 
helium clusters. 

We note that rotation of the rigid body molecule is not included in the path 
integral studies summarized here. While we know that the effect of molecular rota- 
tion on the helium density is limited for slow rotors, this effect is likely to be much 
more important for fast rotors, as was already seen with DMC densities for HCN. 9 
Introduction of the molecular rotation in our path integral computations is there- 
fore planned for future studies. Moreover, the weak coupling angular momentum 
constraint J ~ j appears to be invalid for HCN, at least for the small cluster sizes, 
as demonstrated by the marked inaccuracy of fixed node calculations employing 
free molecule nodal surfaces. All these points, plus the fact that spectra obtained 
by the exact POITSE method using the same projector are nevertheless different 
for the two linear molecules HCN and OCS, underline the marked differences be- 
tween the rotational dynamics of the fast linear rotor HCN from that of the heavier 
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molecules SF 6 and OCS in helium clusters. 
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